Ginzburg-Landau theory of crystalline anisotropy for bcc-liquid interfaces 
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The weak anisotropy of the interfacial free-energy 7 is a crucial parameter influencing dendritic 
crystal growth morphologies in systems with atomically rough solid-liquid interfaces. The physical 
origin and quantitative prediction of this anisotropy are investigated for body-centered-cubic (bcc) 
forming systems using a Ginzburg-Landau theory where the order parameters are the amplitudes of 
density waves corresponding to principal reciprocal lattice vectors. We find that this theory predicts 
the correct sign, 7100 > 7110, and magnitude, (7100 — 7no)/(7ioo + 7110) ~ 1%, of this anisotropy 
in good agreement with the results of MD simulations for Fe. The results show that the directional 
dependence of the rate of spatial decay of solid density waves into the liquid, imposed by the crystal 
structure, is a main determinant of anisotropy. This directional dependence is validated by MD 
computations of density wave profiles for different reciprocal lattice vectors for {110} crystal faces. 
Our results are contrasted with the prediction of the reverse ordering 7100 < 7110 from an earlier 
formulation of Ginzburg-Landau theory [Shih et at, Phys. Rev. A 35, 2611 (1987)]. 
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I. INTRODUCTION 

The advent of microscopic solvability theory 01 Q 0] 
in the 1980s lead to the prediction that the anisotropy 
of the excess free-energy of the crystal-melt interface is 
a crucial parameter that determines the growth rate and 
morphology of dendrites, which shape the microstruc- 
tures of many commercial metallic alloys. This predic- 
tion was largely validated by phase- field simulations 0,0 
of dendritic solidification during the 1990s. More recent 
work in the present decade has focused on the quantita- 
tive prediction of both the magnitude and the anisotropy 
of the interfacial free-energy 7 using molecular dynamics 
(MD) simulations 11 0, 0, H E O E E 13 13 13 
Ej E> Ei E> Ei E| • In parallel, experimental progress 
has been made to determine this anisotropy in metallic 
systems from accurate equilibrium shape measurements 
[23L E> El- which extend pioneering measurements of 
this anisotropy in transparent organic crystals [3 E] • 

MD-based methods, including the cleaving technique 
P Jl l l 13 and the capillary fluctuation method (CFM) 
\WL El , have been successfully developed to compute 7 
and to accurately resolve its notoriously small anisotropy 
of the order of 1%. These methods have been applied to 
a wide range of syste ms, including several elemental face- 
centered-cubic (fee ) ITol ITTl IT3. Hrl E] , body-centered- 
cubic (bcc) [ltilril l2fj | T and hexagonal-close-packed E 
metals, as well as one fee metallic alloy EU modeled with 
interatomic potentials derived from the embedded-atom- 
method (EAML the Lennard- Jones system [3, ^3 1 hard- 
sphere 0, Ej E| and repulsive power-law potentials [j| , 
and, most recently, a bcc molecular organic succinonitrile 
E3 used extensively in experimental studies of crystal 
growth patterns. 

In systems with an underlying cubic symmetry, the 



magnitude of the crystalline anisotropy has been tradi- 
tionally characterized by comparing the values of 7 corre- 
sponding to {100} and {110} crystal faces. MD calcula- 
tions have yielded anisotropy values (7100 —7110)/ (7100 + 
7no) = 0.5 — 2.5%, and experimental values extracted 
from equilibrium shape measurements fall generally 
within this range. What determines physically the pos- 
itive sign and the magnitude of this anisotropy, how- 
ever, remains unclear. One interesting clue is that 
MD-calculated anisotropies generally depend more on 
the crystal structure than on the microscopic details of 
inter-molecular forces for the same crystal structure, and 
anisotropies tend to be consistently smaller for bcc than 
for fee [2(| ■ A striking example of the former is the fact 
that MD studies of Fe pH , and succninonitrile E3 > 
with the same bcc structure but entirely different inter- 
molecular forces, have yielded comparable anisotropies 
around half a percent. The weaker anisotropy of 7 for 
bcc compared to fee is also consistent with experimen- 
tal measurements of anisotropy values in the range of 
0.5 — 0.7% and 2.5 — 5% for the bcc and fee transparent 
organic crystals succinonitrile and pivalic acid, respec- 
tively mil. 

The fact that anisotropy appears to depend more 
strongly on crystal structure than inter-molecular forces 
suggests that it may be possible to predict this critical 
parameter from a continuum density wave description of 
the solid-liquid interface, which naturally incorporates 
anisotropy because of the broken symmetry of the solid. 
The simplest of such descriptions is the Ginzburg-Landau 
(GL) theory of the bcc-liquid interface developed by Shih 
et al. 29]. The order parameters of this theory are the 
amplitudes of density waves corresponding to the set of 
principal reciprocal lattice vectors {Ki}, and the free- 
energy functional is derived from density functional the- 
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ory (DFT) [3(], El with certain simplifying assump- 
tions. 

This theory has yielded predictions of 7 for various bcc 
elements that are in reasonably good agreement with ex- 
periments. Moreover, it has provided an elegant analyti- 
cal derivation of the proportionality between 7 and the la- 
tent heat of fusion, in agreement with the scaling relation 
7n~ 2 / 3 = aL proposed by Turnbull [331 ] and recently cor- 
roborated by MD simulations 0,0], where n is the solid 
atomic density and L is the latent hear per atom. This 
theory, however, predicts the wrong ordering 7 10 o < 7no- 
This apparent failure could be due to the truncation of 
larger \K\ modes, i.e. density waves corresponding to 
reciprocal lattice vectors K with \K\ > \Ki\. However, 
paradoxically, a strong dependence of the anisotropy on 
larger K modes would be hard to reconcile with the weak 
dependence of this quantity on details of inter-molecular 
forces seen in MD studies, since these forces dictate the 
amplitudes of these modes in the crystal where the den- 
sity is sharply peaked around atomic positions. 

To shed light on this paradox, we revisit the simplest 
GL theory of the bcc-liquid interface based on the min- 
imal set of principal reciprocal lattice vectors. Our cal- 
culation differs principally from the one of Shih et al. 
[29j in the derivation of the coefficients of the gradient 
square terms in the GL free-energy functional. Each term 
measures the free-energy cost associated with the spatial 
variation, in the direction n normal to the solid-liquid 
interface, of a subset of equivalent density waves with 
the same magnitude of the direction cosine Ki ■ n and 
hence the same amplitude. Shih et al. choose these co- 
efficients to be proportional to the number of principal 
reciprocal lattice vectors in each subset. This procedure, 
however, turns out to yield an incorrect directional de- 
pendence (i.e. dependence on Ki ■ h) of the rate of decay 
of density waves into the liquid. Here we find that the 
inclusion of the correct directional dependence, as pre- 
scribed by DFT, yields the correct ordering 7100 > 7110 
and a reasonable estimate of the magnitude of anisotropy. 
Furthermore, we validate this directional dependence by 
MD computations of density wave profiles. 

II. GINZBURG-LANDAU THEORY 

We review the derivation of the GL theory and com- 
pare our results to MD simulations in the next section. 
The theory is derived in the DFT framework where the 
free energy of an inhomogeneous system is expressed as 
a functional F — F[n(r)] of its density distribution n(r), 
which can be expanded in the form 



n(r]=n„ l + ^{rl e lR ^ 



(1) 



where the order parameters u^s are the amplitudes of 
density waves corresponding to principle lattice vectors 
(110) of the reciprocal fee lattice, and the contribution of 



larger K mode denoted by "... " is neglected. Expanding 
the free-energy as a power series of the Uj's around its 
liquid value Fi = F[n ] yields 



AF = n ° k B T \ j dra 2 c tj m Uj S j}. + g. 



+b Y. Cl 



(2) 



dui 



dz 



where we have defined AF = F — Fi and the gradient 
square terms arise from the spatial variation of the or- 
der parameters along the direction normal to the inter- 
face parameterized by the coordinate z. The Kronecker 
delta 5 m:n , which equals or 1 for m 7^ n or m = n, 
respectively, enforces that only combinations of princi- 
pal reciprocal lattice vectors that form closed polygons 
Ki + K j + ■ ■ ■ = contribute to the free-energy func- 
tional. Closed triangles generate a non-vanishing cubic 
term that makes the bcc-liquid freezing transition first 
order. The multiplicative factors aj and b are introduced 
since it is convenient to normalize the sums of the c's 
to unity (i.e. J2i c i = h Ef,j ■ c v S o,K t +Kj = etc )- To 
complete the theory, one needs to determine all the coef- 
ficients appearing in the GL free energy functional. 

The coefficients of the quadratic terms are obtained 
from the standard expression for the free-energy func- 
tional that describes small density fluctuations of an in- 
homogeneous liquid 



AF 



drdr* 5n(f) 



-cqf-?\) 



where Sn(f) = n(r) — n and C(\r — r*|) is the direct 
correlation function whose fourier transform 



C{K) = n / dfC(\f\)e 



-iK-r 



(4) 



is related to the structure factor S(K) = [1 — C(K)}^ 1 . 

The two expressions for AF, Eqs. and ©, can 
now be related by assuming that the amplitudes of den- 
sity waves vary slowly across the interface on a scale 
~ 1/ K max where K max is the value of K corresponding 
to the peak of the structure factor. Accordingly, Ui(z') 
can be expanded in a Taylor series about z 



8n(r 



i 

1 d 2 Ui(z) , 



, . duAz) , , , 



2 dz 2 



dz 



(5) 



where the contribution " . . . " involving higher-order 
derivatives can be neglected. Namely, terms propor- 
tional to (z' — z) n d n Ui(z)/dz n ~ l/{K. max w) n , where 
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w is the characteristic width of the solid-liquid interface, 
i.e. the scale over which order parameters vary from a 
constant value in the solid to zero in the liquid. Hence, 
these terms vanish at large n under the assumption that 
w ^> 1/K max . Substituting this expression in Eq. © 
and carrying out the integral over r*, we obtain 



■ Y,\c"{\K i \){k i .- 



dui 


2" 


dz 





(6) 



where C"{K) = d 2 C(K)/dK 2 . Comparing Eqs. © and 
©, we obtain at once 



a 2 Cij 



S(\K 1W \Y 



bc i = -~C"(\K no \)(K i -h) 2 , 



(7) 



(8) 



where we have used the fact that all reciprocal lattice 
vectors have the same magnitude \Ki\ = \Kuq\. Sum- 
ming both sides of Eq. J7J) and using the normalization 



CL2 = 



E 



12 



S(\K 1W \) S(\K 110 \y 



(9) 



and = 1/12. Similarly, summing over i both sides of 
Eq. (jSJl, and using the normalization a = 1, yields 

b=-\Y, C "(\KM)(K l -n) 2 = -2C"(\K lw \), (10) 



where the second equality can be shown to be indepen- 
dent of the direction of h, and 
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c i = -(K i -hf 



(11) 



To make the difference between the above derivation 
of the Cj's and the one of Shih et al. explicit, con- 
sider one of the {110} crystal faces with n pointing 
in the [110] direction. The set of 12 principal recip- 
rocal lattice vectors {Ki} corresponding to the (110) 
directions can be separated into three subsets with 
the same value of (Ki ■ h) 2 : subset I with 8 vec- 
tors ([Oil], [Oil], [Oil], [101], [101], [101], [011], [101]) and 
(K t ■ h) 2 = 1/4, subset II with 2 vectors ([110], [110]) 



and (Ki 



= 1, and subset III with 2 vectors (110, 



[110]) and (Ki ■ n) 2 — 0. Density waves in a given subset 
have the same amplitude denoted here by u, v, and w for 
subsets I, II and III, respectively. 

It follows that the correct coefficient of the gradient 
square terms for a given order parameter u, v, or w, is ob- 
tained using the expression for Ci given by Eq. (jllfl with 



the corresponding value of (Ki ■ h) 2 for the correspond- 
ing subset I, II, or III, respectively. These coefficients are 
Ci = 1/16 for subset I, Ci — 1/4 for subset II, and Ci = 
for subset III. These coefficients yield the gradient square 
terms -C" (\K lw \)\du/dz\ 2 and ~C"(\K lw \)\dv/dz\ 2 in 
the GL free energy functional J2J since there are 8 equiv- 
alent reciprocal lattice vectors in subset I and 2 in sub- 
set II, respectively. The coefficient of \dw/dz\ 2 vanishes 
since principal reciprocal lattice vectors in subset III are 
orthogonal to n and c, = 0. 

In contrast, Shih et al. choose the c^'s to be equal 
for all subsets with a non-vanishing direction cosine 
(subsets I and II), and Ci = for subsets with prin- 
cipal reciprocal lattice vectors orthogonal to n (subset 
III). Since there is a total of 10 reciprocal lattice vec- 
tors in subsets I and II, the normalization condition 
J2i c i — 1 yields Ci = 1/10. These coefficients yield 
the gradient square terms — (8/5)C"(\Ku \)\du/dz\ 2 and 
~(2/5)C"(\K lw \)\dv/dz\ 2 in the GL free energy func- 
tional J2J , which are weighted proportionally to the num- 
ber of reciprocal lattice vectors in each subset, and differ 
from the correct terms derived above. 

For the {100} and {111} crystal faces, the weighting 
procedure of Shih et al. and Eq. (|llf) give coincidentally 
the same coefficients of the gradient square terms. Thus 
these cases need not be repeated here. The results for 
the different crystal faces are summarized in Table [I] 

The determination of all the other coefficients in the 
GL free energy functional is identical to the calculation 
of Shih et al.. The coefficients of the cubic and quar- 
tic terms, cy* and cyy, respectively, are determined by 
the ansatz that all polygons with the same number of 
sides have the same weight, which yields cy^ = 1/8 and 
Cijki = 1/27; for quadratic terms, this ansatz reproduces 
the result cy- = 1/12 derived above since there are twelve 
two-sided polygons formed by the principle reciprocal lat- 
tice vectors. Using these coefficients and identifying each 
Ui with the order parameter u, v, or w, depending on 
whether the corresponding Ki on one side of a polygon 
belongs to subset I, II, or III, respectively, Eq. @ re- 



TABLE I: Comparison of coefficients of square gradient terms 
et predicted by Eq. JTTJ (DFT) and Shih et al. [|<J for the 
{100}, {110}, and {111} crystal faces. For each orientation, 
the 12 principal reciprocal lattice vector are grouped into sub- 
sets where Ki ■ h have the same magnitude in each subset. 
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The corresponding expression of Shih et al. differs by 
the coefficients of \du/dz\ 2 and \dv/dz\ 2 that have an 
extra multiplicative factor of 8/5 and 2/5, respectively, 
as discussed above. Their expressions for AF for the 
{100} and {111} crystal faces are identical to ours since 
Eq. lfTT|) and the equal weight ansatz yield coincidentally 
the same c,'s for these faces. 

Finally, the coefficients (23 and 04 are determined by 
the constraints that the equilibrium state of the solid is 
a minimum of free energy, dAF/du.i\ Ui=U3 — 0, where u s 
is the value of all the order parameters in the solid, and 
that solid and liquid have equal free energy at the melt- 
ing point, AF(u s ) = 0. These two constraints yield the 
relations 0,3 = la^jug and 04 = 0,2 /u 2 s that determine 03 
and (Z4 in terms of a 2 given by Eq. ©, which completes 
the determination of all the coefficients. 



III. RESULTS AND COMPARISON WITH MD 
SIMULATIONS 

The order parameter profiles for {110} crystal faces 
were calculated by minimizing AF given by Eq. i|12|) 
with respect to the order parameters u, v and w, and by 
solving numerically the resulting set of coupled ordinary 
differential equations with the boundary condition u = 
v = w = u s in solid and u = v = w = 0xsi liquid. The 
value of 7110 was computed using Eq. l(T2*|) with these 
profiles. The same procedure was repeated for the {100} 
and {111} crystal faces. 

We used input parameters for the GL theory computed 
directly from the MD simulations in order to make the 
comparison with these simulations as quantitative and 
precise as possible. These parameters include the peak 
of the liquid structure factor ps S'dFuol), which yields 

a 2 = 3.99 using Eq. ©, C"(\K 1W \) = -10.40 A 2 , and 
the amplitude of density waves corresponding to principal 
reciprocal lattice vectors in the solid u s — 0.72. 

The MD simulations were carried out using the EAM 
potential for Fe from Mendelev, Han, Srolovitz, Ackland, 
Sun and Asta (MH(SA) 2 ) [34| and the same thermody- 
namic ensemble and geometries as in Ref. ^J, which 
need not be repeated here. The main difference of the 
present simulations is the way in which the MD results 
were used to calculate density wave profiles. In Ref. 
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FIG. 1: Comparison of planar density profiles n(z) from MD 
simulations (solid line) and the present Ginzburg-Landau the- 
ory (dashed line) for {110} crystal faces. 



the amplitudes were computed by averaging over 
many configurations the instantaneous value of a planar 
structure function (i.e., the magnitude of the complex 
Fourier coefficients of the density). With this approach 
the amplitudes of density waves saturate to a small non- 
vanishing value in the liquid. These amplitudes, however, 
are generally expected to vanish in the liquid that has no 
long range order, consistent with the GL theory. 

To calculate amplitudes that vanish in the liquid, the 
following procedure is followed. We first compute the av- 
erage number density n(r) = n(x, y, z), with z measured 
from a fixed reference plane of atoms in the solid. During 
the MD simulations, only those configurations where the 
solid-liquid interface has the same average position along 
z were considered. As described in detail by Davidchack 
and Laird, |35| this procedure avoids an artifical broad- 
ening of the density profiles due to either the natural 
fluctuations in the average position of the interface or 
Brownian motion of the crystal. The interface position is 
found by first assigning to each atom an order parameter 
proportional to the mean square displacement of atoms 
from their positions on a perfect bec lattice. (The or- 
der parameter calculation is the same as that used in the 
capillary fluctuation method and is described in more de- 
tail in reference [HI). Then an order parameter profile 
as a function of z is computed by averaging within the 
x — y plane and the interface position is that value of z 
where the averaged order parameter is midway between 
the bulk liquid and bulk solid values. 

Next, we compute the x — y averaged density 



n{z) 



L X Ly J 



dxdyn(r), (13) 



which is illustrated in Fig. ^ Lastly, we calculate the 
amplitude of density waves from the fourier transform 



1 



cL x pL y pZj + i 

L x L y Az J Jq J z 



dxdydz n(r) exp(iKi ■ r), 
(14) 
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TABLE II: Comparison of interfacial free energies for different crystal faces (in erg/cm 2 ) and anisotropy parameter £4 = 
(7100 — 7no)/(7ioo + 7110) predicted by Ginzburg-Landau theory with input parameters from MD simulations for Fe with the 
EAM potential of MH(SA) 2 Q (Table ITTTI. and obtained from MD with the MH(SA) 2 potential and two other potentials. 





100 


110 


111 


e 4 (%) 


MD (ABCH) (Ref. [16J) 


207.3 (10.1) 


205.7 (10.0) 


205.0 (10.0) 


0.4(0.4) 


MD (Pair) (Ref. []£]) 


222.5 (14.1) 


220.2 (14.0) 


220.8 (14.0) 


0.5(0.5) 


MD (MH(SA) 2 ) (Ref. [16]) 


177.0 (10.8) 


173.5 (10.6) 


173.4 (10.6) 


1.0(0.6) 


GL (present calculation) 


144.26 


141.35 


137.57 


1.02 


GL (Shih et al. [29]) 


144.26 


145.59 


137.57 


-0.46 



J 



where Zj and Zj+i correspond to sequential minima of 
n(z) and Az = -Zj+i — Zj. In addition, Ui is evaluated at 
the midpoint of this interval, (zj + z J+ i)/2. The order 
parameters u and v were computed for {110} crystal faces 
using Kuo and ifioi, respectively. 

The results of the present GL theory are compared 
to those of Shih et al. and MD simulations in Fig. |3 
and Table [H] Using Eq. (|12|l with the c,'s given by 
Eq. (|ll|l . we obtain the correct ordering of interfacial 
free energies 7100 > 7no an d a weak capillary anisotropy 
(7100 - 7no)/(7ioo + 7110) ~ 1%, consistent with the 
results of MD simulations for bcc elements 0, UjJ 0, 0, 
|20| . while the ansatz of equally weighted q's of Shih et 
al. (with the values listed in Table UJ) gives the reversed 
ordering 7 10 o < 7iio- Note that the predictions of GL 
theory are to be compared to the MD results with the 
MH(SA) 2 potential in TableHTIsince this potential is used 
here to compute input parameters for this theory given 
in Table IIIII MD results for the other potentials are 
mainly included to illustrate the dependence of 7 and its 
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FIG. 3: Comparison of analytically calculated linearized or- 
der parameter profiles u and v for (110) crystal faces near 
the liquid from the present GL theory (solid line) and the GL 
theory of Shih et al. 29] (dashed line), and computed from 
MD simulations using Eq. 1141 with K101 and Kuo for u and 
v, respectively (solid circles). 




100 110 120 130 140 150 



Z(A) 

FIG. 2: Comparison of numerically calculated nonlinear or- 
der parameter profiles u and v for (110) crystal faces obtained 
from the present GL theory (solid line) and the GL theory of 
Shih et al. |29| (dashed line) and computed form MD sim- 
ulations using Eq. 1141 with K\oi and Kuo for u and v, 
respectively (solid circles). 



anisotropy on details of interatomic forces. 

Fig. n shows that the planar density profile predicted 
by GL theory, obtained by substituting Eq. with 
the numerically calculated order parameter profiles into 
Eq. I|13|) . is in remarkably good agreement with MD 
simulations on the liquid side. The discrepancy on the 
solid side is due to the fact that GL theory neglects the 
contribution of larger \K\ reciprocal lattice vectors that 
contribute to the localization of density peaks around bcc 
lattice positions in MD simulation. 

Fig. [21 shows that the amplitude profiles and the in- 
terface widths predicted by GL theory are in good agree- 
ment with MD simulations. The MD results clearly val- 
idate the directional dependence of the rate of spatial 
decay of density waves in the liquid that is the main de- 
terminant of the anisotropy of 7. This directional depen- 
dence is most clearly seen by examining the amplitude 
profiles on the liquid side of the interface. In this region, 
the amplitudes of density waves are sufficiently small that 
one can neglect the cubic and quartic terms in the GL 
free energy functional. The resulting linear second order 
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differential equations for u and v obtained by minimizing 
this functional can be solved analytically, and have ex- 
ponentially decaying solutions that are compared to the 
MD results in Fig. |3| The coefficients of the gradient 
square terms in the free energy functional control the de- 
cay rates. The u and v profiles (i.e. the amplitude of 
density waves corresponding to K\oi and Kuo) calcu- 
lated with coefficients that depend on the angle between 
the principal reciprocal lattice vectors and the interface 
normal through Eq. which is consistent with DFT, 

have different decay rates that are in good quantitative 
agreement with the MD results. In contrast, u and v 
profiles calculated based on the ansatz of equal weights 
for the Cj's [2t| have the same spatial decay rate, which 
does not agree with the MD results. 

It is interesting to note that Mikheev and Chernov 
(MC) j2EH3> m a formulation of the anisotropy of the 
solid-liquid interface mobility, also stress the importance 
of the decay rate of the amplitude of density waves. The 
MC model predicts crystal growth rates and anisotropics 
that are in qualitative agreement with MD simulations 
of FCC systems. The theory, however, is linear in the 
sense that only the effective widths of the density profiles, 
which arc allowed to vary with K and n, are required and 
the authors make no attempt to compute, as was done 
here, the full amplitude profile as a function of z. 

Finally, even though we have focused primarily in this 
paper on crystalline anisotropy, it is useful to re-examine 
the prediction of GL theory for the magnitude of 7 and 
for the Turnbull coefficient using input parameters from 
the present MD simulations. Shih et al. [2^| derived 
an analytical expression for the magnitude of 7 in the 
isotropic approximation where all the order parameters 
are assumed to have the same profile through the inter- 
face, i.e. Ui(z) = u for all i. In this approximation, 
the free energy density reduces to the sum of the gradi- 
ent square term b\du/dz\ 2 and a quartic polynomial in u. 
The stationary profile u(z) that minimizes the free en- 
ergy is then an exact hyperbolic tangent profile and the 
analytical expression for the interfacial energy is 



n k B T„ 
6 



(15) 



Furthermore, Shih et al. related the latent heat (per 
atom) to the temperature variation of the inverse of the 
peak of structure factor proportional to a 2 (Eq. |5J , 
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where N is the number of atoms in the system. This 
yields the expression for the Turbull coefficient 
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(a 2 b) 



1/2 



L 3T m dd2/dT\ T=1 



(17) 



which Shih et al. evaluate using parameters for the hard 
sphere system [2||. Using the values of the various co- 
efficients obtained from MD simulations listed in table 



IIIII Eq. (|15fl yields a value of 7 = 147.4 erg/cm 2 in rea- 
sonably good agreement with the average values of 7 for 
the different crystal faces in Table ITD obtained from MD 
simulations and the fully anisotropic GL calculation with 
different order parameter profiles. With the same coef- 
ficients, Eq. Ijl6(l yields a latent heat value L = 0.114 
eV/atom about 30% lower than the MD value in Table 
IIIII where the difference can be attributed to the con- 
tribution of larger K modes that are neglected in GL 
theory. Eq. l|17fl i n turn predicts a value of the Turnbull 
coefficient a — 0.45 that is about 25% larger than the 
MD value a m 0.36, owing to the underestimation of the 
latent heat of melting in GL theory with input param- 
eters of Table IIIII from the present MD simulations. In 
the future, it would be interesting to test how the Turbull 
coefficient predicted by GL theory (Eq. ITTjl varies with 
input parameters computed from MD simulations using 
different interatomic potentials. 



IV. CONCLUSIONS 

We have revisited the simplest GL theory of the bec- 
liquid interface whose order parameters are the ampli- 
tudes of density waves corresponding to principle recip- 
rocal lattice vectors. We find that, despite its simplicity, 
this theory is able to predict the density wave structure 
of the interface and the anisotropy of the interfacial en- 
ergy, in reasonably good quantitative agreement with the 
results of MD simulations. 

A main determinant of the anisotropy of the interfacial 
energy in this theory is the rate of spatial decay of den- 
sity waves in the liquid. This decay rate must depend on 
the angle between principal reciprocal lattice vectors and 
the direction normal to the interface for this theory to be 
consistent with DFT. This directional dependence, which 
we validated quantitatively by MD simulations, is a di- 
rect reflection of the underlying crystal structure. There- 
fore, the present results provide a simple physical picture 
of the strong relationship between crystal structure and 
crystalline anisotropy, consistent with the findings of a 
growing body of MD-based and experimental studies of 
crystalline anisotropy. 

An interesting future prospect is to extend the GL 
theory to other crystal structures, and in particular fec- 
liquid interfaces. This requires, however, to consider the 
coupling of density waves corresponding to the principal 
reciprocal lattice vectors to larger K modes, which makes 
the theory intrinsically more complicated. 
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da 2 /dT (K" 1 ) 


u s 7 (erg/cm 2 ) 


L (eV/atom) 


MD (MH(SA) 2 ) 
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0.00163 


0.72 175(11) 


0.162 
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